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Abstract  -  In  this  study,  we  present  a  method  to  perform  a 
numerical  simulation  of  the  flow  dynamics  of  the  Cerebrospinal 
Fluid  (CSF)  based  on  anatomical  Magnetic  Resonance  Images 
(MRI).  The  Computational  Fluid  Dynamics  (CFD)  software, 
written  in  language  C,  integrates  different  numerical  schemes  to 
solve  the  governing  equations.  The  time  derivatives  were 
discretized  using  the  Crank-Nicolson  scheme.  The  equation  of 
continuity  was  modified  by  introducing  an  artificial 
compressibility  and  discretized  by  a  finite  difference  scheme.  The 
meshed  boundary  of  the  CSF  was  immersed  in  a  Marker-And- 
Cell  staggered  grid  for  to  take  into  account  fluid-structures 
interactions.  Equations  of  hydrodynamics  were  solved  with  an 
iterative  method  under  different  quasi-static  conditions.  The 
anatomical  basis  of  our  simulations  was  generated  from 
individual  MRI  scans.  The  surface  of  the  anatomical  flow 
channels  of  interest  was  extracted  by  segmentation  and 
triangulated.  In  parallel  to  the  acquisition  of  the  anatomical  data 
CSF  flow  has  been  measured  by  MRI.  To  characterize  a  whole 
cardiac  cycle  sixteen  equidistant  velocity  measurements  have 
been  performed.  In  addition,  a  home  made  software  was 
implemented  to  visualize  computed  data  (velocities,  pressure). 
Keywords  -  CSF  flow;  Computational  Fluid  Dynamics;  Magnetic 
Resonance  Imaging;  Intra-cranial  dynamics;  CSF-Brain 
interactions. 

I.  Introduction 

The  aqueduct  of  Sylvius  is  the  narrow  channel  connecting 
the  third  and  fourth  ventricle  (Fig  1).  It  is  a  narrow  pathway 
for  the  Cerebrospinal  Fluid  (CSF)  flow  and  hence  may  be  a 
privileged  site  for  development  of  hydrocephalus. 

The  geometry  of  the  aqueduct  of  Sylvius  has  often  been 
studied  [1,2]  and  numerous  models  of  the  physics  of 
hydrocephalus  have  been  developed  [3].  Jacobson  [4]  realized 
a  numerical  simulation  of  the  CSF  with  a  commercial 
Computational  Fluid  Dynamics  (CFD)  software  that  couldn’t 
take  into  account  elastic  and  compressive  properties  of  the 
surrounding  brain  tissue. 


Fig.  1  :  three  dimensional  representation  of  intra-cerebral 
CSF  spaces. 


Fig.  2  :  3D  model  of  the  aqueduct  of  Sylvius. 

The  meshing  is  made  of 458  vertices  and  881  triangles. 

The  aim  of  this  study  is  to  realize  a  realistic  numerical 
simulation  of  the  CSF  flow  including  fluid- structure 
interactions. 

II.  Methodology 

A.  Anatomy  modeling 

The  first  part  of  this  work  comprises  the  generation  of  a 
numerical  model  of  the  brain-CSF  interface  (called 
parenchyma)  in  which  the  flow  is  simulated.  Three 
dimensional  T2-weighted  MR  Images  were  acquired  on  a  1.5 
T  Signa  scanner  (GEMS,  Milwaukee)  to  provide  a  good 
contrast  between  brain  structures  and  CSF.  A  software  was 
realized  using  IDL  5.4  (Interactive  Data  Language,  RSInc 
France)  in  order  to  process  image  data.  The  desired  CSF 
spaces  were  isolated  from  surrounding  tissues  in  three 
dimensions  by  virtual  cutting  of  the  upper  and  lower  end  and 
subsequent  application  of  a  region  growing  based 
segmentation  algorithm.  This  structure  was  meshed  by 
triangulation  and  the  resulting  meshing  was  organized  to  fully 
describe  the  relationships  between  triangles,  vertices  and 
edges. 

B.  Numerical  simulation 

The  CFD  code  was  written  in  ANSI  C  language  based  on 
the  work  of  M.  Fluck  [5].  Under  normal  physiological 
conditions,  the  CSF  behaves  as  an  incompressible  Newtonian 
fluid  in  a  laminar  flow  [6].  Density  and  viscosity  were  p  =  103 
kg.m"3  and  p  =  10"1  Pa.s  respectively  [4]. 

Under  these  conditions,  the  governing  laws  of 
hydrodynamics  are  the  equation  of  continuity  : 

V  .0  =  0  (1) 

and  the  Navier-Stokes  equation  : 
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p —  +  p(uS/)u  =  -Vp  +  //At]  +  f 
dt  (2) 
where  u  represents  the  velocity  vector  (whose  components 
are  u,  v  and  w  for  the  x,  y  and  z  directions,  respectively),  p 

the  pressure  and  f  the  forces  that  act  on  the  fluid. 

Different  numerical  schemes  have  been  used  to  solve  these 
equations.  For  the  linearization  of  temporal  terms  of  the 
Navier-Stokes  equations,  the  Crank-Nicolson  scheme  was 
chosen  for  to  provide  both  a  good  convergence  and  acceptable 
truncation  errors.  The  Artificial  Compressibility  Method  was 
applied  with  the  equation  of  continuity  and  then  linearized 
with  a  backward  difference  scheme. 

This  system  of  equations  was  solved  for  a  finite  number  of 
points.  The  discretization  of  the  spatial  domain  was  realized 
with  the  multi-grid  Immersed  Boundary  Method.  For  this 
method  [7],  the  numerical  model  of  the  parenchyma  was 
immersed  in  a  three  dimensional  Cartesian  grid.  Grid  points 
were  labeled  ‘Inner’  or  ‘Outer’  depending  on  their  location 
inside  or  outside  the  aqueduct.  Then,  the  flow  was  calculated 
for  inner  points.  For  each  point  of  the  grid,  velocity 
components  and  pressure  were  staggered  (Fig.  3)  following 
the  Marker-And-Cell  (MAC)  method. 

After  the  linearization  step  of  relations  (1)  and  (2),  an 
iterative  method  was  used  to  compute  a  numerical  solution  for 
each  grid  point  from  the  following  algebraic  equation  system 
whose  unknowns  were  velocity  u(UjVjW)  and  pressure  p  : 

Ui[n+'1-UiP+>(.Lu,(up,\^1,\V{n,p[n)=0  (3) 

pf+'-pf  +  LD  (u jn,v[n,wf1)=0 
where  m  represents  the  iteration  index  and  /  the  position  of  the 
point  in  the  MAC  grid.  In  this  expression,  the  linearized 
equations  are  written  as  functions  :  Lui  stands  for  Navier- 
Stokes  for  the  component  number  i  (z  =  1,  2,  3)  of  velocity; 
and  D  stands  for  the  divergence  operator  in  the  equation  of 
continuity.  X  and  k  are  two  parameters  that  have  to  satisfy  the 
following  conditions  to  provide  the  convergence  of  the 
iterative  method  : 


0  <  k  <  (Ax)2- 


0  <  A< 


(Ax)2 

2k 


2  (Ax)2 

Re+  2  At 

1  (Ax)2 

Re+  4At 


(4) 


Re  is  here  the  Reynolds  number,  At  the  time  step  and  Ax  the 
grid  spacing.. 

A  correct  solution  was  obtained  for  iteration  index  m+1  if : 

\ur'-ur\<  £  < <  i  (5) 

for  each  of  the  variables  u,  v,  w  and  p. 

Following  the  Immersed  Boundary  Method  proposed  by 
Peskin  [7],  the  coupling  between  flow  and  membrane 
deformation  was  insured  by  a  three  dimensional  Dirac  delta 
function.  Membrane  forces  were  calculated  with  the  Finite 
Element  Method  on  each  vertex  of  the  meshed  surface.  The 
membrane  had  to  mimic  the  influence  of  the  brain  tissues 
surrounding  the  aqueduct.  Consequently,  an  approximation 
has  been  made  and  the  mechanical  properties  taken  for  the 


Fig.  3  :  2D  staggered  arrangement  of  velocity 
components  and  pressure. 


TABLE  1 

Mechanical  constants  for  different  material. 


Material 

Young’s 
Modulus  E 
(N.ni-2) 

Poisson  ratio  v 

reference 

Gray  Matter 

50 

0.49 

[81 

White  Matter 

5 

0.49 

[81 

“homogeneous 

Brain” 

10 

0.35 

[3] 

Density 

(kg.m-3) 

Viscosity 

(Pa.s) 

CSF 

103 

10'1 

[41 

%  of  the  Cardiac  Cycle 

Fig.  4  :  distribution  of  the  CSF  velocities  over  a  Cardiac 
Cycle  (CC). 

Positive  values  mean  downward  (caudal)  direction  while 
negative  mean  upward  (cranial)  direction. 


membrane  were  values  the  of  Young’s  Modulus  E  and 
Poisson  ratio  v  found  in  the  literature  [3,8]  for  the  brain 
supposed  homogeneous  and  visco-elastic  (cf.  table  1). 
However,  the  software  imposed  the  membrane  to  be  rigid  in  a 
first  validation  step. 

C.  Quasi-static  flow 

The  simulation  was  performed  under  quasi-static  conditions 
over  a  cardiac  cycle,  i.e.  the  input  of  the  simulations  have 
been  the  velocities  acquired  by  measurements  of  CSF  flow 
made  by  Phase  Contrast  MR  images  over  a  cardiac  cycle  [9]. 
“Input”  velocities  are  the  boundary  conditions  imposed  on  top 
of  the  model.  Sixteen  simulations  have  been  run.  For  each 
run,  the  input  velocity  profile  was  chosen  to  be  uniform  to 
reflect  the  constriction  between  the  3rd  ventricle  and  the 
aqueduct. 
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Fig.  5  :  Velocity  norm  (iso-contour)  and  vector 
field  (black  lines)  along  the  aqueduct  at  about 
80%  of  the  Cardiac  Cycle. 

Further  software  based  on  IDL  has  been  developped  to 
analyse  and  visualise  results  of  the  numerical  simulation  for 
both  flow  velocities  and  pressures. 

III.  Results  and  Discussion 

Fig.  1  and  2  show  the  segmented  structures  and  the  meshing 
resulting  from  the  image  processing  step,  respectively.  This 
step  allows  the  representation  of  thin  anatomical  structures 
even  with  complex  geometry.  The  quality  of  the  model 
depends  on  image  resolution.  With  a  resolution  of  about  0.7 
mm,  a  sufficient  geometry  was  computed. 

As  the  Fig.  5  shows,  velocities  develop  along  the  aqueduct 
in  a  parabolic  profile.  Computed  pressures  (Fig.  6)  are  in  good 
agreement  with  CSF  dynamics  [9],  i.e.  at  80%  of  the  Cardiac 
Cycle,  lower  pressure  in  the  3rd  ventricle  allow  the  “fill  flow” 
to  develop. 

Different  approximations  have  been  made  concerning  the 
brain  homogeneity,  isotropy  and  visco-elasticity.  These 


Fig.  6  :  Pressure  distribution  along  the  aqueduct  at 
about  80%  of  the  Cardiac  Cycle. 


hypotheses  are  rather  coarse,  but  led  to  a  simpler 
implementation  of  the  program.  No  flux  through  the 
membrane  has  been  allowed.  A  porosity  term  in  the  model  of 
the  brain,  as  proposed  by  Tada  et  al.  [8],  was  omitted. 

IV.  Conclusion 

The  presented  method  allows  to  simulate  interactions 
between  brain  and  CSF  flow  based  on  anatomical  data. 

This  study  is  the  first  step  toward  a  hilly  dynamic 
simulation  with  moving  boundaries.  Further  implementation 
will  include  calculation  of  strain  in  the  membrane  and, 
therefore,  allow  a  better  description  of  the  intra-cranial 
physics. 

Acquisition  of  CSF  velocities  in  hydrocephalic  patients  and 
individual  modeling  and  simulation  with  these  data  may 
improve  our  understanding  of  the  pathogenic  mechanisms  of 
hydrocephalus  development.  Moreover,  influence  of 
geometry  (e.g.  partial  or  total  aqueduct  obstruction)  or  brain 
mechanical  properties  on  hydrocephalus  may  be  another  point 
of  interest. 


References 

[1]  R.  J.  Last  and  D.  H.  Tompsett,  ’’Casts  of  the  cerebral 
ventricles,”  J Br  Surg,  vol.  164,  pp.  525-563,  1953. 

[2]  S.  S.  Kemp,  R.  A.  Zimmerman,  L.  T.  Bilaniuk,  D.  B. 
Hackney,  H.  I.  Goldberg,  and  R.  I.  Grossman, 
’’Magnetic  resonance  imaging  of  the  cerebral  aqueduct,” 
Neuroradiology ,  vol.  29,  pp.  430-6,  1987. 

[3]  G.  Tenti,  J.  M.  Drake,  and  S.  Sivaloganathan,  ’’Brain 

biomechanics:  mathematical  modeling  of 

hydrocephalus,”  Neurol  Res ,  vol.  22,  pp.  19-24,  2000. 

[4]  E.  E.  Jacobson,  D.  F.  Fletcher,  M.  K.  Morgan,  and  I.  H. 
Johnston,  ’’Fluid  dynamics  of  the  cerebral  aqueduct,” 
Pediatr  Neurosurg,  vol.  24,  pp.  229-236,  1996. 

[5]  M.  Fluck,  ’’Dreidimensionale  modellbildung  und 
simulation  von  roten  blutzellen  im  hydrodynamischen 
scherfeld,”  in  Diplomarbeit  in  Physic.  Aachen:  RWTH, 
1998,  pp.  138. 

[6]  A.  K.  Ommaya,  ’’Mechanical  properties  of  the  nervous 
system.,”  J Biomech,  vol.  l,pp.  127-138,  1968. 

[7]  C.  S.  Peskin,  ’’Numerical  analysis  of  blood  flow  in  the 
heart.,”  J  Comput  Phys,  vol.  25,  pp.  220-252,  1977. 

[8]  Y.  Tada  and  T.  Nagashima,  ’’Modeling  and  Simulation 
of  Brain  Lesions  by  the  Finite-Element  Method.,”  IEEE 
Engineering  in  Medicine  and  Biology ,  pp.  497-503, 
1994. 

[9]  O.  Baledent,  M.-C.  Henry-Feugeas,  and  I.  Idy-Peretti, 
’’Cerebrospinal  Fluid  Dynamics  and  relationship  with 
blood  flow  :  an  MR  study  with  semi-automated  CSF 
segmentation,”  Investigative  Radiology ,  vol.  July  -  in 
press,  2001. 


